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BODIES 

by  Zhu  Youlan,  Wang  Ruquan,  and  Zhong  Xichang 

Chinese  Academy  of  Sciences,  Institute  of  Computer  Technology 


1.  INTRODUCTION 


From  the  latter  part  of  the  50’ s  to  the  present,  many  methods 
have  appeared  for  finding  a  numerical  solution  to  the  problem  of 
supersonic  inviseid  flow  about  blunt-nosed  bodies.  Of  these,  one  type 
is  classified  as  the  non-stationary  methods  and  the  other  type  as  the 
steady-state  methods.  Among  the  steady-state  methods  there  is  also 
the  finite  difference  method,  the  integral  relationship  method,  the 
straight-line  method,  etc.  All  of  these  methods  relate  to  smooth 
objects  and  all  can  obtain  satisfactory  results.  But,  since  the  non- 
stationary  method  must  undergo  a  process  of  stabilization  with  respect 
to  time,  the  machine  time  expended  is  very  great.  In  regard  to  the 
finite  difference  method,  in  order  to  obtain  sufficiently  accurate 
results,  it  is  required  that  grid  points  cannot  be  too  few,  and  thus 
machine  storage  capacity  must  be  rather  great.  Also  a  large  amount 
of  computer  time  is  expended.  In  comparison,  the  straight-line 
method  has  many  advantages.  For  example,  its  algorithms  are  simple 
and  the  required  storage  capacity  is  small,  and  using  very  few  rays 
makes  it  possible  to  obtain  satisfactory  results.  This  paper  mainly 
describes  our  work  in  the  field  of  supersonic  flow  about  blunt-nosed 
bodies  using  the  straight-line  method  of  calculation. 
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We  used  the  straight-line  method  to  conduct  extensive  calculations 

of  supersonic  flow  about  blunt-nosed  bodies.  The  calculated  objects' 

shapes  were  ellipsoids  with  various  axial  ratios  as  well  as  objects 

ancormrig 

which  were  nearly  disc-shaped.  The  Mach  numbers  of  theAflow  ranged 
from  1.5  to  «•.  Under  conditions  of  axial  symmetry,  calculations 
were  carried  out  on  not  only  frozen  gas  with  y  =  1.4  but  also  on  flow 
in  equilibrium  air  and  non-equilibrium  air. 

Furthermore,  we  also  used  the  straight-line  method  to  calculate 
flow  in  the  supersonic  region  and  flow  with  a  sharp  angle  of  attack. 

We  conducted  many  checks  along  the  way  on  the  calculated  results. 
All  of  the  checks  indicated  that  the  results  obtained  using  the 
straight-line  method  were  quite  satisfactory. 

2.  FORMULATION  OF  THE  PROBLEM 


2.1  The  basic  equation 

Inviscid  and  non-heat-conducting  steady-state  gaseous  flows  were 
taken  into  consideration.  The  aerodynamic  equation  in  the  spherical 
coordinate  system  (r,  0,  is 


~i\  (pur*  6)  +  —  (p»»r  sin  6 )  +  ~~  (p*or)l  — 

r  sin  6  'Or  00  0q>  i 


dtt  r*  +  tv*  .  1  dp 

- 1 - -r1-  “  0 

it  r  p  Qr 

*L  +  !2L-2i°wt  +  JL*Z.-o 

it  r  r  Tp  06 

^  +  "Z  +  51livw  +  _J_  it  -  0 

dt  r  r  pr  s.n  $  dtp 

~  -  0 

it  it 


(2.1) 


where  t-  “  “  ~~~  +  *  —  +  *• — ~ ,  respectively,  are  the  velo- 

dt  dr  r  do  r  ba  9  Otp 

city  components  in  the  directions  of  r,  0,  tp  ,  p  is  pressure,  p  is 
density,  and  c  is  sonic  velocity.  All  the  quantities  in  the  equation 
are  dimensionless  quantities  and  their  dimensional  factors  are, 
respectively , 

p  —  pmVi,  P  —  Pm,  V-,  r  ~  R, 
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Here,  subscript  00  represents  the  quantity  of  oncoming  flow  and  R0  is 
the  radius  of  curvature  of  the  apex  of  the  object’s  nose.  In  order  to 
simplify  calculation,  we  introduce  the  following  coordination  conver¬ 
sions  : 


f 


_ f  -  c(fl>  ?>  e 

no.v)-c(e,  9) 


0.  9 
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where  r  —  G(0,  9)  and  r  —  F{6,  9)  are  the  equations  for  the  object  surface 
and  shock  wave,  respectively.  It  is  evident  that  in  the  ({,0,9) 
coordinate  system  the  shock  wave  lies  within  the  £  =  1  plane  and  the 
object  surface  lies  within  the  5  =  0  plane.  Let 

*  -  -  <c#  +  (et),  fi~-  -4-  (c,  +  fO.  «  —  F  —  C 

OD  6 

and  since  JL  —  JLJL  0  «  d  a 

dr  "  *  d{  *  a«  "  «  9{  as 

_L.  jL  -  i.  ,  1  a 

tin  a  69  e  S'  tin  0  69 

consequently,  equation  (2.1)  can  be  written  as 


,$lL  +  p(r%!L  +  0§!L  +  fi&!:)- Ftt  *  +.L|L_F, 

a*  \  a$  df  /  as  p  df 

•  —  4-  —  —  F,,  .  ^  -  F.,  f, 

Si  p  a*  a«  p  a*  vas.  ajV 


or  written  in  the  solved  form  of  ....  £2!  : 

0f  *  a,c 

■  * 

dp  r_  ^[F,*  —  p(.tF,  4-  oF,  -4-  FF«)1  *♦*  Ft* 


a? 


■*  —  r*el 


a*  f>l  «  a*  J 

~  —  —  f  F, - -  ^-1 

aj  -  I  p  a*  J 

—  —  —  [  f,  —  — 

af  *  i  e  es  j 

es  •  l  p  a{  J 


(2.2) 


(2.3) 


where 


*  "  +  po  +  wfi 

r*  -  r‘  +  «»  +  p* 
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2.2  BOUNDARY  CONDITIONS 


(1)  Shock  wave  conditions 

Based  on  the  shock  wave,  they  must  satisfy  the  shock  wave  relat¬ 
ionship 


(2.4) 


where  quantities  with  subscript  oo  represent  wave-front  quantities, 
quantities  without  subscripts  represent  wave-tail  quantities,  Y*  re¬ 
presents  the  projection  of  velocity  in  the  direction  of  the  shock  wave 
normal,  h  is  enthalpy,  and  (nt,  n2 ,  n3)  are  the  respective  directional 
cosines  of  the  shockwave  normal,  i.e., 
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(2)  Object  surface  conditions 

At  the  object  surface  they  must  satisfy  normal  velocity  for  zero 
conditions ,  i . e . , 

f  *  u c  —  »C(  —  w  - *  0  (2.5) 

tin  6 

3.  NUMERICAL  SOLUTION 

In  order  to  carry  out  numerical  solution,  we  introduo*-^  a  number 
of  rays  in  the  solution  region.  For  example,  in  the  6  direction  we 
introduced  some  coordinate  planes  with  0  =  9^  =  const,  in  the  q  direc¬ 
tion  we  introduced  some  coordinate  planes  with  <p  =  =  const  and  took 

the  intersection  lines  of  these  coordinate  planes  as  the  rays.  For 
fixed  Ci  we  used  the  value  of  the  flow  parameter  based  on  the  ray  as 
the  nodal  point  value  to  form  an  interpolation  polynomial  and  thus 
were  able  to  correspondingly  determine  the  partial  derivatives  for  $  and 

4.  Letting  equation  (2.3)  be  established  on  the  basis  of  the  rays 
will  lead  to  a  set  of  ordinary  differential  equations  and  the  corre¬ 
sponding  boundary  value  problems  will  become  boundary  value  problems 
of  a  set  of  ordinary  differential  equations.  For  this  we  converted 
the  boundary  value  problems  into  initial  value  problems  and  found  the 
solution  by  iteration,  i.e.,  given  the  shockwave  form,  we  determined 
the  flow  parameters  of  the  wave  tail  from  shockwave  conditions  (2.4) 
and  using  these  as  the  initial  values,  integrated  the  set  of  ordinary 
differential  equations  and  then  checked  whether  or  not  the  object  sur¬ 
face  flow  parameters  thus  obtained  satisfied  object  surface  conditions 
(2.5).  If  not,  we  then  readjusted  the  shockwave  form  until  the  object 
surface  conditions  were  satisfied.  The  usual  methods,  such  as  the 
four-step  Runge-Kutta  method,  can  be  used  for  integrating  the  set  of 
ordinary  differential  equations.  Some  of  the  actual  treatments  are 
described  below. 

3.1  Equations  based  on  the  9=0  axis 

Assuming  a  flow  field  with  ^  =  0,  the  IT  plane  is  symmetrical  and 
the  0=0  axis  is  totally  within  the  symmetrical  plane.  From  equation 


5 


(2.2)  it  can  be  seen  that  since  sin  ^appears  in  the  denominator, 
equations  with  0=0  must  also  be  treated.  Let  *•■»,»({,  o,  0),  and  it  is 
evident  that  *({,  o,  -  *•({)«>» Vt  o,  ~  9.  *•({.  0,  *)-«({,  0,  0), 

p({,  o,  ?)— p({,  o,  o),  p(f,  o,  <?)-  p((,  o,  o).  We  used  Lobida's  law  for  the  per¬ 
centage  terms  appearing  in  equation  (2.2)  and  were  easily  able  to  ob¬ 
tain  an  equation  based  on  0  =  0.  In  principle,  it  is  well  to  select 
equations  based  on  a  plane,  but  due  to  errors  in  numerical  calcula¬ 
tion,  for  different  the  results  obtained  can  be  different.  To  eli¬ 
minate  such  differences  we  integrated  these  equations  for  <p  from  0  to 
IT  and  obtained  the  desired  equations.  For  example,  after  multiplying 
the  third  equation  in  (2.2)  by  cos  and  subtracting  the  fourth 
equation  multiplied  by  sin  <?,  we  then  carried  out  integration,  and 
since 


f*  /  dr  du>  .  \  , 

l  t  I - cot  ®  —  — —  »n  <?  )  a  q 

J*  \bi  6f  J 


(  (r*  ■+■  at  fit v")  —  dq> 

J*  d;  . 


6vm 


Of 


1  .  S§Pj  ■*  dp 

J  («ce,  -  «,  tin 

|  (F,cotq>  —  F,  tin  —  —  ~{j^  [»'*  “  tin  ycotgjj 

+  —  cot  «pl  4<p  +  —  Mr*l  F* 

f  ae  J  2  i 

consequently,  we  can  obtain 

-•  +  ®L  -  f? 

Of  p  Of 


* 


similarly,  we  can  obtain 


5«.  + JLS£.-F? 

d{  p  di 


■•t  +  '(f+”’  If) " Fr 


written  in  the  form  of  •  •  • ,  ,  it  is  easy  to  obtain  the  desired 

df  Of 

equations 


?JL- 

fMFN*  - 

Of 

®£.  - 

ILK 

Off 

l  «* 

6u  ^ 

J.  [  f  J  —  . 

Off  " 

.•r* 

Os! 

Of 

.•r* 

,*t  —  T»v 


V  (2.6) 


where 


.•-ij  UtQitpJty 
«%  -  r«  o*^ 


T**  -  r>  +  0 


*> 


f r  - 

f?  - 

f:  - 

f?  - 


-  t !"'  j!f«  “***'  +  'S.’lj '* + »"’! 

“7 1'"  i.’ I?  “*»'»- 7  -“} 

2 1  f  f  ■  f  ,  /  5t»  .  8w  .  \  1  5»  1 

- I  *  l  TT  eo*‘T  “  T  *m  <p  cot  f  cot  u>  4u> 

*  U*l  \06  fi9  /  p  dti  J 

(•(i!L-.fla£\co  ^ 

*  J*  \6d  SO  / 


*  *1 

-  Mt<  ? 

2  1 


3.2  Interpolation  polynomials 

Since  it  is  assumed  that  the  flow  field  for  ^  =  0  and  the  x  plane 
are  symmetrical,  therefore  a  solution  is  required  only  between  0  £  «f 
^tf.  We  introduced  k  +  1  half-planes  between  0  and  Tt.  At  the  same 
time,  n  +  1  cones  are  formed,  0  =  §K  -  const  (0#  =  0),  A=  C,  1,  . . . ,  n.  They 
intersect  half-planes  4  =  4  *,  and  ^  s  efi  +Tf  forming  (2n  +  1)  rays.  Note 
symmetry  of  flow  parameters  and  that  those  for  ^  =  are  equal  to 

those  for  *  =  +  ft ,  or  differ  in  sign.  Consequently,  for  a  fixed  5, 

the  values  for  2n  +  1  rays  at  4  =  <?i  and  4  =  tr  -  can  be  used  to  form 
0  interpolation  polynomials  of  the  2n-th  degree: 

(2.7) 

<-0 

’  r 

g  represents  a  flow  parameter.  To  save  time  i:;  calculation,  we  did 
not  precalculate  coefficient  a’t ,  but  used  the  following  algorithm: 
since  the  interpolating  function  and  its  derivative  can  both  be  re¬ 
presented  as  a  linear  combination  of  functional  values  for  the  inter¬ 
polation  nodal  point,  and  since  U  j  linear  combination  coefficient  is 
related  to  both  the  position  of  the  interpolation  nodal  point  and  the 
position  of  the  interpolation  point,  therefore,  once  these  nodal 
points  and  interpolation  points  are  fixed,  the  coefficients  can  be 
determined  and  when  the  function  and  derivative  values  of  the  various 

•T 

points  must  be  calculated,  it  is  only  necessary  that  the  functional 
values  for  these  coefficients  and  nodal  points  be  scalar  products. 
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Taking  equation  (2.7)  as  an  example,  assume  that  the  nodal  point  in 
0m(m  =0,  1,  ...»  2n).  We  must  calculate  g(8)  at  some  interpolation 
point  0.  Since  at  the  nodal  point  there  is  the  condition 

22  •  “  (.  (m  -•  0,  1,  •••,  2h) 


gm  ■■  /($„).  Or  written  in  the  form  of  a  matrix: 


Ma  =  g 


where 


i  fi,  6>- 

1  o, . 


i  eim 


Therefore,  we  have 

a  =  iHg 

Thus,  the  equation  for  g(0)  can  be  considered  as 

g(0)  =  d*  •  a 

rewritten  as 

*<*)-d’  •  (M-g) 

-  (M—  d)'  ■  t 

b*  .  g  (2.8) 

where  b  =  M“1#d,  and  d*  =  (1,0,  .  ..,02  n).  it  is  evident  that  when 

the  nodal  points  and  interpolation  points  are  fixed,  M*  and  d  can  both 

be  determined  and  thus  we  will  also  be  able  to  obtain  b.  Likewise, 

since  .Li, 

x#(0)  —  22  ••'(’i)'”' 

(•I 

—  dr  • » 

-<w—  «*.)•» 

where  d?  =  (0,  1,  20,  ...,  2n02n-1),  the  prime  (')  represents  a  deri¬ 
vative  of  9  and  therefore  we  consider  the  derivative  can  be  similarly 
treated . 
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With  respect  to  0,  when  0  is  fixed,  we  use  the  value  where  it  in¬ 
tersects  the  k  +  1  half-plane  to  construct  the  trigonometric  interpol¬ 
ation  polynomial  for  0.  The  method  mentioned  ubove  is  also  used  in 
calculation.  For  even  functions  we  take 


therefore 

we  now  have 


COlV 


M 


1  1  1 
„  |  cos  9.  cc«  qp,  cos  9* 


cos^9,  eos^9,»**  co»*9^ 


•*  (0,  —  sin  9,  *ia 9 co* *“'9) 

**  “*  (/*.  Xt»  Ik) 


for  -i- 1  tJq>,  —  I  gcot a>d<p  we  can  also  write  similar  expressions  whi 

*  Ji  tr  Ji 

have  net  been  enumerated.  For  odd  functions  we  take 


ch 


hence 


Now 


t-1 

X  ■“  53  «, cos*  9  sin  9 
<  ■  0 


d,)*f 


*/•  - 


s»n  9, 
tin  9,  cos  9, 


>«n  9,cos*“,9i 


sin  9,- 


sin  9)eos^“'9» 


sin  9»-i 


sin  94-,cos<",94-/ 


di,*  “  (cos  9,  —  sin*9  4-  eof*9»  •••,  —  2)c»s,*",9*>n,9  4*  cos<-19) 

In  addition  to  the  method  described  above,  we  also  used  the  fol¬ 
lowing  method  to  construct  interpolation  polynomials.  For  even  func¬ 
tions  of  we  take  .  , 

c(i)  -  53  23  x.-(*)®‘«o*,9 

•  •  Q  #  »  0 


for  odd  functions  we  take 

«■({)  —  (53  23  «'.l(*)0ieo«v)  sin  9 

'  • •  ©  I  • 0  ' 

In  order  to  bring  the  function  and  its  derivative  for  0  to  a  single' 
value  at  0  -  0,  we  must  add  the  appropriate  conditions  to  the  poly- 
nomial.  Now  to  explain  such  conditions. 


# 


9 


For  example,  for  a  shock  wave  form  of  *  F($,  9),  when  0  s  0,  we 

naturally  require  that  it  be  independent  of  $.  This  merely  requires 
that  =  0  when  j  4^  0.  At  the  same  time,  in  order  to  ensure  that 
the  normal  has  a  definite  direction,  when  0s 0  the  direction  of  the 
normal  under  spherical  coordinates  should  have  a  form  such  as  this: 

(«,  *to»9>,  —  Aun 9),  where  a  and  b  are  arbitrary  constants.  This  requires 
that  Fi,  —  0  (j  1)  .  Therefore,  for  F  the  polynomial  should  be 

«  < 

K®.  9)  •  +22  F.t&cot'tp  (2.9) 

•  »  I  »  •• 


It  is  evident  that  p,  p  3hould  also  take  the  same  form  as  F.  The 
above  conclusion  and  the  v  and  w  forms  can  still  be  obtained  by  the 
following  means.  We  let 

g  -  <«<{)  +  SZ 

<•1  , 

■  k 

1  •  1  i  «e 

*>  —  —  9  +  (E  E  V 

where  g  can  represent  p,p,u,F,  and  the  characteristics  »>({,  0,  <p)  —  ^coiy, 
“•({»  0,  sO  “  —  ^*«n v,  were  used,  and  therefore  we  have  i(N-o(,  1), 

^  —  0(/  **  0),  —  «*■«.  Entering  the  above  formula  into  equation  (2.2), 

notice  that 

—  c(0,  q> )  —  G#(0,  C)co*<p, 

00 

-£2.1  «—  —  C*(0,  0)*in  <p 

iinfllw 


(where  we  assume  that  the  object  surface  is  at  =  0  and  the  if  plane 
is  symmetrical.  Therefore,  when  Q 0,  we  can  obtain 


+  'cP.^f  -  P.1  C;  +  {(F„  -  G.))^  -  Jim  F, 

,/jS?  +  Slits  -  lim  F, 

Poe  il 

(st  4*  -  G.~  +  ZiiSZG*}  i£i)  COi,r  _  lim  F, 

-  U  ^  -  c<  +  i(r"  ~  Cft)  sin  v  -  lim  F. 

\  pm  dl  >  »-• 


lim  F, 

(-• 


lim  r , 

l-l 


10 


where 


•t  "  —  ««l  ct.  +  £(f„  —  <7*,)] 

'i  -  C(c,  0)  +  £lf»  -  G( 0,  0)] 

fi- 

Me 

The  right  side  of  the  above  formula  is  written  as  a  polynomial  of 

C03  <f>  (or  multiplied  by  sin  &  and  using  the  linear  independence  of 

1 ,  cosy,  •  •  •  ,co*V'  it  can  be  easily  deduced  that  g  should  take  the  form  of 

equation  (2.9)  and  v  and  v  should  take  the  form 

«  < 

f  *“  4-  t >t!6  •+■  o^q>  4-  ^ i\.8  cotfip 

f-i  l  • u 

/  .*  \ 

x  «■»  —  i',,  tin  9-  —  vtjO  cos  q>  tin  <p  4-  (  «’(l0'cos'<p  J  tin  <p 

',«»  i«0  ' 

At  the  same  time,  we  have  also  obtained  an  equation  based  on  0  =  0: 


®C  «»  tff 


“  (F«  ~  G(0,  0))l t/wp„  4-  pw(k,;  4-  2i/l#  4-  2«i*)j 


d  M, 


p* 

!*£«<  _  c*  +  ‘(Fn  ~ 

Pot 


( Fgj  C(0>0))rol(»l,  t'sx) 


el 


“  (Fg,  —  G(0,  0))  [iti(t',0  4-  t>„  +  ug,)  4-  &1] 

1  Pt«J 

••  -  -  C*.  -  c( o.  o))  (,„  -  u*  „,) 


(2.10) 


3-3  Iteration  methods 

Just  as  mentioned  in  the  beginning  of  this  section,  we  solved  the 
boundary  value  problems  by  iteration.  In  fact,  this  can  also  be 
understood  by  solving  the  problem  of  a  set  of  nonlinear  transcendental 
equations,  i.e.,  selecting  an  Fj  set  (representing  the  shock  wave  form 
at  the j-th  ray),  we  use 

•  •  •,  F„)  0  (/  •»  1,  2,  •••  j  m) 

where  qx  =  0  ,  i.e.,  boundary  condition  (2.5).  We  used  the  Newton 
method,  i.e.,  the  change  in  F|  value,  6Fj  ,  satisfies  the  following 
equation : 
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V  M.  eF,  - 

,  r  °  r  ' 

•  •  0  ^  * 


—  q,  0  —  l.  •••.  ">) 


Normally,  the  partial  derivative  is  incapable  of  being  expressed 

in  analytical  form,  therefore  we  used  the  numerical  method,  i.e., 


fV  _  q  ( Ft,  •  •  • ,  f F ,  +  r.-,,  •  •  • ,  —  SXF,,  •  •  • , 

df,  AF, 

Using  this  method,  each  time  we  iterate,  we  then  must  integrate  m  +  1 

times,  hence  an  extraordinary  expenditure  of  machine  time.  To  save 

time,  the  simplified  Newton  method  can  be  used,  but  since  is 

oF. 

probably  not  very  exact,  the  convergence  rate  may  be  relatively  slow. 

To  make  ~~  a  little  more  accurate,  but  without  increasing  the  cal- 
oF, 

culated  value  very  much,  we  suggest  the  following  method.  Let  Q  re¬ 
present  a  vector  consisting  of  qt . qm,  and  let  R  represent  a 

vector  consisting  of  .  . .  ,  Fm.  The  Q  to  R  guidepost  [sic]  was  express¬ 
ed  as  Q’.  Given  that  Q(R0)  =  Q0  and  m  points  R4,  ...,  Rm  are  located 
outside  the  vicinity  of  R0  and  it  is  known  that  Q(Ri)  =  Qi ( i  =  1,  ...» 
m),  and  if  R0-Ri(i  =  1,  2,  ...»  m)  is  linearly  independent,  then  QJ 
can  be  approximated  from  Rj  and  Qi .  In  fact,  since 

P.  Pt  +  £>,( R,  —  A,)  («  —  1 ,  2,  •  •  • ,  m) 


F-)  -S.(F„ 


we  therefore  have 


(Pi  P«.  pa  P()  pl(  K I  R»,  •  •  •  .  F  m  A») 


thus  we  have 


P»  ***  ( P>  ~  P*»  *  * '  >  P-  “  P*X *•»  ‘  *  * » R*)  * 
and  therefore  can  obtain  the  iteration  formula 


—  A'.  —  (F._.  —  Rm,  A,_,  —  R.XQw-m  —  P.»  Q.-i  ~~  P.X'P. 

( 

*  «■  m  +  1 ,  •  •  • 


(2.11) 


It  is  obvious  that  using  the  above  formula  to  carry  out  iteration, 
first  to  find  Qt,  ...,  Qm+l  from  Rj,  ...,  Rra+1  requires  integration 
m  +  1  times,  but  subsequent  iteration  only  requires  integration  one 
time.  Now  a  simple  estimate  ^or  convergence  rate.  It  is  obvious 


that  when 


-(A,  -  •  E 
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equation  (2.11)  is  a  differencing  Newton  formula  in  which  E  is  the 
unit  matrix  and  is  a  scalar  quantity.  S*nce  we  now  have 

Qm* I  (  pi  '  P<B*|»  ■  *  t  Qm  '  P»«|)  “  +  o(AF*) 

therefore  we  have 

p;%',  -  ^ _ .Oi’,,  +  -  Hm*tQm*,  +  «(Af) 

consequently  for  the  differencing  Newton  formula  we  have  the  estimat¬ 
ing  formula 

R •  -  R..,  -  R*  -  Rm.,  +  . 

“  R*  —  Rm~t  +  Pf.~JlPi.4l  +  (^«4lP«»l  Pfn.'OP— *» 

-  «(1!R*  -  R„*,II;)  +  *<  a  Flip-,  I!) 

or  written  as 

HR*  -  RJ  <  A ||R*  -  R..,B»  4-  BAfilp._,H  («-  1,  2,  •  ••) 
where  R*  represents  the  genuine  solution,  j|  |(  represents  a  modulus, 
and  A  and  B  are  both  appropriate  constants.  Likewise,  for  a  simpli¬ 
fied  differencing  Newton  method  there  ir,  the  estimating  formula 


!iR*  -  R.t!  <  W  -  +  CAF|!p..,||  +  CRB.*,  -  R.-,H lip.-, D 

(»■*»•♦*  2,  ••  *)  " 

For  the  method  we  proposed  there  is  the  estimating  formula 


i !**  --*.51  <  •*!; Rm  -  R.-,!i  +  f)!lp:l,D  .uP  HR..,  -  R.-.-.li'Hp.-itl 

,<•  <• 


where  A,  B,  C,  and  D  are  all  constants. 


It  is  not  difficult  to  see  that  when  &f  is  on  the  same  order  or 
smaller  than  f|  R*  -  Rn_l  ([  ,  the  differencing  Newton  method  basically 
maintains  the  convergence  rate  of  the  Newton  method.  But,  generally 
II  -  «n-1 11  increases  as  n  increases,  and  therefore  the  converg¬ 
ence  rate  of  the  simplified  Newton  method  is  relatively  slow.  Since 
p;L,!!  mp  "R.-,  —  R.-,-.li  —  oO),  and  !UP  liR.-,  —  R. decrease  with  an  in- 
crease  in  n,  therefore  the  method  we  propose  can  converge  somewhat 
faster  than  the  simplified  Newton  method.  Actual  calculations  have 
confirmed  this  point. 


3.4.  Initial  shock  wave  selection  and  object  surface  value  extrapol¬ 
ation 

When  using  the  method  mentioned  above  to  conduct  iteration, 
whether  the  shock  wave  is  well  selected  or  not  is  rather  important  for 


13 


calculating  whether  or  not  it  will  be  successfully  carried  out.  For 
this  reason,  we  used  the  results  already  obtained  on  the  basis  of 
step-by-step  interim  measures  of  certain  parameters.  For  example, 
when  we  want  to  consider  the  results  of  flow  under  a  certain  angle  of 
attack  n#»  angle  of  attack  n  can  be  selected  as  a  parameter.  After 
the  results  of  nc  have  been  calculated  (for  example  n0  =  O  »  then  we 
can  regard  it  as  the  initial  value  of  n0  +  An.  and  after  calculating 
the  results  of  n0  +  An.  we  can  then  use  the  results  of  n0 »  n0  +  An 
and  use  linear  interpolation  to  determine  the  initial  value  of  n0  + 
2An-  And  so  forth  and  so  on,  until  the  results  of  n*  have  been  deter¬ 
mined.  As  to  first  calculation  of  the  required  initial  shock  wave 
form,  this  can  be  accomplished  by  using  available  results. 

It  can  be  seen  from  equation  (2.3)  that  a  =  0  on  the  object  sur¬ 
face  and  therefore  integration  is  not  carried  out  up  to  the  object 

surface.  In  order  to  determine  the  object  surface  value,  we  used  the 
extrapolation  method,  i.e.,  integration  is  carried  out  until  a  ce^tairn 
£*  is  reached  (for  example  £»  =  0.1),  then  using  the  £  value  of  sever¬ 
al  adjacent  points  we  extrapolate  the  object  surface  value,  for  exam¬ 

ple,  using  the  values  £  =  0.3.  0.2,  0.1  is  regarded  as  two-step  inter¬ 
polation.  It  should  be  mentioned  in  passing  that  if  another  form  of 
calculation  is  used,  for  example  the  implicit  form,  to  carry  out  in¬ 
tegration,  extrapolation  can  be  completely  avoided.  For  conditions  of 
axial  symmetry,  we  designed  a  mode  and  carried  out  the  calculations 
and  the  results  indicate  that  it  was  a  success.  It  is  not  described 
in  any  detail  here. 

4.  RESULTS  OF  CALCULATIONS 

We  used  the  methods  described  above  to  conduct  extensive  calcul¬ 
ations  on  axisyrametric  and  three-dimensional  flow  about  blunt-nosed 
bodies.  The  shapes  of  the  calculated  objects  were  ellipsoids  and 
nearly  disc-shaped  objects  with  various  axial  ratios  (the  objects  are 
represented  by  the  equation  z"  +  {x*  +  y*)*S**  /*  ,  and  n^  2).  The  Mach 
number  range  of  the  oncoming  flow  was  1.5  ^  o©  .  For 
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conditions  of  axial  symmetry,  in  addition  to  frozen  gas  with  y  =  1.4, 
we  also  conducted  calculations  for  equilibrium  air  and  non-equilibrium 
air.  The  method  described  in  [12]  was  used  to  calculate  the  thermo¬ 
dynamic  functions  of  equilibrium  air.  See  our  other  articles  for  non¬ 
equilibrium  gases.  The  accuracy  of  the  computational  results  was 
checked  in  various  ways.  One  was  the  carrying  out  of  the  computations 
themselves,  for  example,  the  use  of  different  numbers  of  rays  and 
different  integral  step  lengths  and  conducting  checks  on  certain 
relationships  which  the  flow  field  must  satisfy,  such  as  the  condition 
of  total  conservation  of  energy.  Another  was  making  comparisons  of 
results  obtained  from  experiments  and  from  other  methods  such  as  the 
integral  relationship  method.  All  the  checks  indicated  that  the 
accuracy  of  the  computed  results  was  quite  satisfactory.  Part  of  the 
computed  results  is  given  below. 


Figure  1  shows  the  shape  of  shock  waves  and  sonic  speed  lines  for 
the  flow  of  frozen  gases  about  a  sphere  at  different  M«o  numbers.  It 
is  evident  from  the  figure  that  the  sonic  speed  line  forms  appear  in 
two  types:  At  M  >  3,  the  limiting  characteristic  line  is  the  second 
family  of  characteristic  lines  to  reach  the  object  surface  sonic 
speed  point;  and  at  M<  3?  the  limiting  characteristic  line  is  made 
up  of  the  first  family  of  characteristic  lines  in  contact  with  the 
sonic  speed  line  (generated  by  the  object  surface)  and  the  second 
family  of  characteristic  lines  (generated  by  the  shock  wave). 


Fig.  1.  The  shock  wave  and  sonic 
speed  line  shapes  for  the  flow  of 
frozen  gases  about  a  sphere 


1.1  1.4  u  1.0 
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Figure  2  shows  the  object  surface  pressure  distribution  for  frozen 
gas  flow  about  a  sphere. 


§ 

io- 


Fig.  2.  Distribution  of  pressures  over 
the  object  surface 


Figure  3  shows  shock  wave  shape  and  sonic  speed  line  location  for 
different  objects,  Moo  =  3,  T  =  I-1*.  For  a  very  blunt  body,  for  ex¬ 
ample  n  }  20 ,  the  position  of  the  shock  wave  will  generally  remain 
unchanged.  With  a  decrease  in  the  radius  of  curvature  of  the  object 
surface  in  the  vicinity  of  the  sonic  speed  points,  for  a  fixed  M  oo 
number,  beginning  with  a  certain  curvature  the  sonic  speed  line  will 
exhibit  an  inflection  point. 


Fig.  3  Shock  wave  locations  and  sonic 
speed  line  shapes  (Mo©  =  3) 


Figure  4  present  equilibrium  gas  flow  about  a  sphere  and  the  shock 
wave  locations  and  sonic  speed  line  shapes  at  different  Mo©  numbers. 
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Pig.  4.  Shock  wave  and  sonic  speed  line 
shapes  for  the  flow  of  equilibrium  gas  about 
a  sphere 

KEY:  a)  equilibrium;  b)  frozen 


The  flow  conditions  are:  Moo  =  4,  poo  =  0.87  X  103  dynes/cm^,  fo©  s 
0.95  X  10'®  g/cm3  and  Moo  a  20,  poo  =  0.122  X  10®  dynes/cm^,  $  oo  = 
0.192  X  10~3  g/cm3.  It  is  evident  from  the  figure  that  dissociation 
causes  the  shock  wave  location  and  sonic  speed  line  to  markedly  differ 
from  that  of  frozen  gas  with  Y  =  1.4.  The  shock  wave  draws  very  near 
to  the  object  surface. 


Fig.  5.  Shock  wave  and  sonic  speed  line 
shapes  for  the  flow  of  non-equilibium  gas 
about  a  sphere 


Figure  5  shows  shock  wave  location  and  sonic  speed  line  sha  for 
the  flow  of  non-equilibrium  gas  about  a  sphere.  The  conditions  of  the 
oncoming  flow  are:  p  ••  =  0.947  X  1 0 3  dynes/cm^,  =  0.123  X  10"5 
g/cm3 ,  R  =  5  cm.  It  should  be  noted  that  the  M  «o  =  20  sonic  speed 
line  has  a  peculiar  shape  at  the  shock  wave. 
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Fig.  6.  Stagnation  point  shock 
wave  detachment  distance  along 
with  changes  in  Mo* 

KEY:  a)  frozen;  b)  equilibrium 
c)  non-equi librium;  d)  dynes/ 
cm2;  e)  g/cm3. 


Figure  6  shows  the  relationship  between  the  Moo  number  and  the 
detachment  distance  of  the  stagnation  point.  For  frozen  flow,  with 
M  >  10,  it  basically  remains  unchanged.  For  equilibrium  air,  with 
an  increase  in  M  #*  ,  a  change  in  detachment  distance  occurs  without 
monotonicity. 

When  calculating  axisymmetric  flows  normally  five  rays  are  used. 
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Fig.  7.  Shock  wave  shapes  and  sonic 
speed  lines‘within  a  symmetrical  plane 
along  with  cnanges  in  angle  of  attack 

KEY:  a)  sonic  speed  line;  b)  shock 
wave;  c)  direction  of  oncoming  flow; 

d)  stagnation  point  (Mo*  =  4,  3); 

e)  Results  from  [12] 
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Figure  7  presents  the  flow  about  an  ellipsoid  with  6  =  1.5  and  the  j 
shock  waves  and  sonic  speed  lines  within  a  symmetrical  plane  at  Hn  : 

3  and  4  along  with  changes  in  angle  of  attack  n*  Figure  9  presents 
the  object  surface  pressure  within  a  symmetrical  plane  along  with 
changes  in  angle  of  attack.  Figure  7  also  presents  the  location  of 
the  stagnation  points  and,  within  the  accuracy  of  the  figure,  the 
stagnation  points  for  M  oo  =  3  and  M  po  =  4  duplicate  each  other.  .  | 

When  the  angle  of  attack  changes,  the  stagnation  point  moves  over  the-  ] 

object  surface  at  a  nearly  uniform  rate.  Figure  8  presents  the  ; 

shock  waves  and  sonic  speed  lines  within  the  ^  =  rt/2  plane  along  with  .j 

changes  in  angle  of  attack.  In  the  figure  it  can  be  seen  that  changes 
in  shock  wave  shape  are  slow,  but  the  changes  in  sonic  speed  lines  are 
relatively  quick. 
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Fig.  8.  Shock  wave  shapes  and  sonic 
speed  lines  with  changes  in  angle  of 
attack  =4,6=  1.5) 
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Fig.  9.  Object  surface 
pressure  within  a  sym¬ 
metrical  plane  along 
with  changes  in  angle  of 
attack  (6  =  1.5) 

KEY:  a)  Results  from  [  12] 
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Fig.  10.  Shock  wave  shapes  and  sonic 
speed  lines  on  a  symmetrical  surface 
(M*o  =  5.8,  n  =  20' 


Figure  10  presents  the  flow  about  a  disc-shaped  object  for  M  ••  = 
5.8  and  n  =  20  and  the  shock  waves  shapes  and  sonic  speed  line  locat¬ 
ions  on  a  symmetrical  plane  at  different  angles  of  attack.  Figure  11 
presents  the  object  surface  pressure  distribution  within  an  axial 
[sic]  plane. 


Fig.  11.  Object 
ure  distribution 
metrical  plane 
(Mm  -  5.8,  n 


surface  press- 
within  a  sym- 

=  20) 


When  calculating  three-dimensional  flow,  we  generally  take  four 
f  planes,  and  within  each  ^  plane  we  take  four  rays.  And  because 
the  axes  are  common,  altogether  we  take  13  rays.  The  coordinate 
system  z-axis  extends  into  the  symmetrical  plane  of  the  flow  field 
and  for  convenience  of  calculation  we  usually  use  the  included  angle 
between  the  flow  field  symmetrical  plane  and  the  body  symmetrical 
axis  (angle  (£  =  (0.5  -  0.6)  >f  ,  and  yj  is  the  angle  of  attack). 
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In  regard  to  the  question  of  whether  or  not  maximum  entropy  value 
is  reached  on  the  object  surface,  it  can  be  seen  from  our  calculations 
that  for  certain  types  of  objects,  within  a  given  range  of  error, 
maximum  entropy  is  reached. 

We  employed  several  means  to  check  the  results  of  our  calculations 
We  used  different  numbers  of  rays,  for  example,  for  axisymmetrical 
flows  we  conducted  separate  calculations  using  five  rays  and  three 
rays  and  the  results  indicated  an  error  of  within  1J.  We  also  employ¬ 
ed  different  integral  step  lengths.  For  example,  for  flow  about  a 
sphere  with  M  oo  ~  6  and  f  =  1.4,  from  the  shock  wave  to  the  object 
surface,  integrals  of  10  steps,  20  steps,  and  80  steps.  The  relative 
error  in  the  results  obtained  with  10  and  20  steps  did  not  exceed 
0.3%  and  between  20  to  80  steps,  there  was  at  least  the  same  three 
significant  digits.  This  also  explains  that  we  did  not  have  to  be 
overly  concerned  about  rounding  off  error  increases. 

We  made  a  comparison  between  the  results  for  M  «o  *  3,  £  =  1.5, 
rj  =  15°  and  the  results  obtained  by  Telenin  [16]  and  it  can  be  seen 
from  Figure  1  that  they  are  exactly  the  same.  For  frozen  flows  about 
ellipsoids  and  spheres  with  $  =  1.5,  the  results  of  our  calculations 
agree  with  those  obtained  by  Belotserkovskiy  using  the  integral  rela¬ 
tionship  method  up  to  three  significant  digits. 

We  also  conducted  integral  checks,  i.e*,  we  also  checked  the  acc¬ 
uracy  of  integral  equation 

jj^pu  •  do  -  o 

JJ-u  •  Jo  +  pXi  •  do  —  0 
(•*-1,2,3) 

Jj#PP"Vu  •  da  —  0 

where  x*  represents  the  unit  vector  in  the  3  axial  directions  under 
rectangular  coordinates,  u^  is  the  projection  of  velocity  loss  u  in 
these  three  directions,  is  any  curved  surface  which  does  not  in¬ 
clude  the  shock  wave.  See  Table  1  for  the  results  of  integration  for 


Mo»  =4,  6  1.5,  n  =  1 0° ,  20°.  Besides  the  total  integral,  Table 

1  also  gives  the  integral  for  the  shock  wave,  for  the  object  surface 
and  for  the  cone.  It  can  be  seen  from  Table  1  that  the  calculated 
results  are  correct  and  that  integrals  for  the  shock  wave  and  cone 
differ  only  in  the  third  digit.  Object  surface  conditions  are  also 
rather  well  satisfied,  being  on  the  order  of  10**1*. 
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»j>fitttt<!© 

0.09011x2 

-0.01979X2 

-0.07HSX2 

0.00033  x2 

20* 

yr>  fittfi60 

0 

0 

0 

0 

« SfittaCj) 

0.092*1  x2 

—  0.30*17x2 

-0.34423X2 

-0.00139x2 

aOO 

-0.056379x2 

0.00001x2 

0.03663x2 

0.00033X2 

Table  1.  Integral  checks  M«e  =  4,  6  =1.5,  f}  =  0.6n 

KEY:  (a)  equation;  (b)  shock  wave  surface;  (c)  object  surface; 

(d)  cone  surface;  (e)  Total;  (f)  mass;  (g)  momentum  in  direction  x 
ion;  (h)  momentum  in  direction  y;  (j)  momentum  in  direction  z; 

(k)  entropy 


Furthermore,  we  calculated  the  total  energy  at  all  the  nodal 
points  and  the  entropy  of  each  nodal  point  on  the  object  surface. 

From  Table  2  it  can  be  seen  that  with  M««  =  V;  <$  -  when  n  ^  15°, 

the  relative  error  of  total  energy  is  less'than  1$ .  From  Table  3  it 
can  be  seen  that  with  M  *•  =  4 ,  6  =  1.5,  and  n  ^  15°,  the  object  sur¬ 

face  entropy  differs  only  by  1  at  the  third  significant  digit. 

Table  2.  Total  energy  difference 
The  exact  value  for  j.  _JC_  JP.  is  0.65625,  at  M*o  =  4 

*>  r-i 


1 

•Oe.K  Cxifi* 

*0 

o-  ; 

0.0010 

0.13% 

3* 

0.0019 

0.30% 

10* 

0.0034 

0.51% 

13* 

0.0031 

0.78% 

20* 

O.00S7 

1.33% 

KEY:  (a)  total  energy  maximum  dev¬ 
iation;  (b)  relative  error 


Table  3-  Entropy  value  of  object  surface 
nodal  points  (:n'-r,"0 

When  M  *•  =  4 ,  6=  1.5,  8  =  0 . 6  n  »  n  =  0 , 

the  exact  value  is  -2.31904 


-2.3U5J 


-2.1IVM 


-2.il/2J 


-2.31517  -2.1(611 

-2.11116  -2.10*4$ 
-2,11121  -2.11*42 
-2.11*11  -2.11141 
-2.J2091  -2.11*11 

-2.J2O10  —  2.31551 
-2.11*12  -2.12117 
-2.31919  j  -2.1171V 
-2.j:o:7  -2.12111 

-2.J1276  -2.102*1 
—  2.  HSU  -2.11141 
-2.11*51  -2.11901 
-2.12U15  -2.32121 


In  short,  using  the  straight-line  method  to  calculate  flow  about 
smooth  bodies,  the  accuracy  of  the  results  is  satisfactory. 

Comrade  Rui  Weiming  participated  in  part  of  the  work,  Comrade 
He  Jiaomin  gave  us  much  assistance,  and  Comrade  Feng  Kang  enthusiast 
ically  directed  our  work  and  for  this  we  express  our  thanks  to  them. 
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